#######################################################################################
# File: MESMLM-ARp.R                                                             
# Author: Xun PANG                                                             
# Creation Date: 4/5/2014                                                            
# Description: 	MCMC simulation for estimating binary MESMLM-AR(p) model
#               based on the research paper, "Varying Responses to Common Shocks and
#               Complex Cross-Sectional Dependence: Dynamic Multilevel Modeling with
#               Multifactor Error Structures for Time-Series Cross-Sectional Data "
#
# Usage:          MulFac.AR <- function(y,X1, Wi, Ai, Ut, TP, Unit.Index, ...)
#                 MulFac.AR0 <- function(y,X1, Wi, Ai, Ut, TP, Unit.Index, ...)  
# Arguments:
#         y     response variable, dichotomous (0/1)
#        X1     matrix of covariates with fixed effects
#        Wi     matrix of covariates with subject-varying effects
#        Ut     matrix of covariates with time-varying effects
#        Ai     matrix of covariates explaining the subject-varying effects
#        TP     number of unobserved factors
# Unit.Index    index of subjects
# Time.Index    index of time periods
# timeint.add   add an time-specific intercept? "TRUE" or "FALSE". Default is "FALSE"
# unitint.add   add an unit-specific intercept? "TRUE" or "FALSE". Default is "FALSE"
#         m     number of iterations to be returned
#    burnin     number of initial iterations which are discarded
#    priors     beta0--- mean vector of the multivariate normal distribution of beta,
#                        all the fixed effects (beta in the reduced form of the model)
#               B0--- variance matrix of the multivariate normal distribution of beta,
#                        all the fixed effects (beta in the reduced form of the model)
#               d0, D0 --- degree of freedom and scale matrix of the Wishart prior
#                          on b_i, the subject-level residual
#               e0, E0 --- degree of freedom and scale matrix of the Wishart prior
#                          on c_t, the time-level residual
#     nlag      the lag order p, integer >= 1
# initial values
#               init.phi --- vector of initial values of the autoregressive coef.
#               init.beta, init.b, and init.c --- initial values of the fixed and
#                   random coefficients
#               init.ystar --- initial values of the latent response variable
# tracking      every "tacking" iterations, return the information about how many
#               iterations in total have been done
#  monitor      a string contains the names of parameters whose MCMC output are
#               will be returned. The string has to be a subset of
#               ("rho", "beta", "bi", "ct", "D", "E", "ystar", "u") which is the default
#     thin      thin the chain by recording each "thin" iterations
#                                       
######################################################################################
  # Call a file with functions necessary for MCMC updating
  source("hidden.function.R")

  # The main function for MCMC updating
  MulFac.AR <- function(y,X1, Wi, Ai, Ut, TP, Unit.Index, Time.Index,  m=10000, burnin=5000, init.phi,
                        init.ystar=0, init.beta=0,unitint.add=0,timeint.add=0, init.b, init.c,
                        init.theta=0, init.gamma=0, beta0, B0, D0, d0, E0, e0,  nlag, 
                        monitor=c("rho", "beta", "bi", "ct", "D", "E", "ystar", "u", "theta", "gamma"),
                        thin=1, tracking=50){
    
     ############################
     #y \in {0, 1} error checking
     ############################
     yy <- y[!is.na(y)]
     if (sum(yy!=0 & yy!=1) > 0) {
       cat("Elements of y equal to something other than 0 or 1.\n")
       stop("Check data and call MulFac.AR() again. \n") 
     }
     
    #****************** Arrange Data for MCMC updating  **********************#
     
     # To use the RegressionData for this model specification, the arguments of
     #  Ft and St are set to be certain values since the data inpute does not contain
     #  information for Ft and St
     Ft <- "NULL"
     St <- as.matrix(rep(1, length(y)))
     # Call RegressionData() and rearrange the data to fit the reduced model
     NewData <- RegressionData(yob=y, Unit=Unit.Index, Time=Time.Index, Fixed=X1,
                               UnitRandom=Wi, TimeRandom=St,
                               UnitPred=Ai, TimePred=Ut)

    Y <- NewData[[3]]
    entry1 <- 1: length(Y)
    newid1 <- NewData[[1]]
    newid2 <- NewData[[2]]
    id11 <- Index(newid1)     #the unit diex
    uniqunit <- as.numeric(table(id11))   # how many time periods for each unit
    jj <- which(uniqunit<= nlag)
    if (length(jj)>0){						
      entry2 <- entry1[id11==jj]
      Y <- Y[-entry2]
      X <- as.matrix(NewData[[4]][-entry2,])
      WP <- as.matrix(NewData[[5]][-entry2,])
    }else{
      X <- as.matrix(NewData[[4]])
      WP <- as.matrix(NewData[[5]])
    }
    WPP <- ncol(WP)
    KW <- ncol(Ut)
    FB <- ncol(X)
    W <- as.matrix(cbind(WP, X[, (FB-KW+1):FB]))
    if (unitint.add==TRUE){
      W <- cbind(1, W)
    }
   N <- length(Y)
     
   # Check the dimension of the parameter vectors
     
   #beta: fixed parameter vector 
   PFB <- length(beta0)
   if (FB!=PFB){
      cat("Length of the prior on fixed effects is incorrect.\n")
      cat("Length of fixed effect should be ", FB, "\n", sep=" ")
      stop("Respecify beta0 and B0, call GLMMARp.Binary() again. \n") 
    }
   # bi: random fixed parameter vector 
   RB <- ncol(W)
   if(is.matrix(D0)){ 
     PRB <- ncol(D0)-TP
   }else{
     PRB <- length(D0)-TP
   }
   if (RB!=PRB){
       cat("Dimension of the prior on unit random effects is incorrect.\n")
       cat("Length of D0 should be ", RB+TP, "*", RB+TP,  "\n", sep=" ")
       stop("Respecify D0, call GLMMARp.Binary() again. \n") 
    }

   # in case of unbalanced TSCS data structure,
   # create new indices to accomendate such a situation 						
   entry <- seq(1: N)                   # the global length of observations
   if (length(jj)>0){
     # unit id repeating t times t=# observations
     id1 <- Index(id11[-entry2])                
     # year id repeating n times n=# unit in each time period
     id2 <- (newid2[-entry2])-min(newid2[-entry2])+1        
     GU <- max(id1)                       # how many units
     GT <- max(id2)                       # how many time periods in total
     id3 <- id1
     GU <- length(unique(id3))
     n=1
     while (n <= length(unique(id3))){
       rows <- entry[which(id3==n)]                   
       nT <- length(rows)
       TT <- id2[rows]
       if (nT!=1){
         for (j in 2:nT){
          GHK <- TT[j]-TT[j-1]
          if (GHK >1){
            PL <- rows[j]
            id3[PL:N] <- id3[PL:N]+1
          }
        }
       }
       n <- n+1
     }

    entry3 <- 1: length(Y)
    uniqunit <- as.numeric(table(id3))
    jj <- which(uniqunit<= nlag)
    entry4 <- entry1[id3==jj]
    if(length(entry4)>0){
      Y <- Y[-entry4]
      X <- as.matrix(X[-entry4,])
      W <- as.matrix(W[-entry4,])
      id1 <- Index(id1[-entry4])
      id2 <- Index(id2[-entry4])
      id3 <- Index(id3[-entry4])
    }
   }else{
     id1 <- id11
     id2 <- Index(newid2)
   }
    id1 <- Index(id1)                 # unit id repeating t times t=# observations
    id2 <- (id2)-min(id2)+1        # year id repeating n times n=# unit in each time period
    id3 <- id1						
    GU <- max(id1)                       # how many units
    GT <- max(id2)
    GUU <- length(unique(id3))


     #############################
     # Storage and Initial Values
     #############################
   total.return <- m+burnin
   if ("rho"%in% monitor){ 
     if (nlag==1){
       phi <- rep(NA, total.return)
       phi[1] <- init.phi
     }else{
       phi <- matrix(NA, nrow=total.return, ncol=nlag)
       phi[1,] <- matrix(init.phi, nrow=1, ncol=nlag)      
     }
   }else{
     if (nlag==1){
       phi <- init.phi
     }else{
       phi <- matrix(init.phi, nrow=1, ncol=nlag)
     }
   }
   
   if ("beta"%in% monitor){ 
     beta <- matrix(NA, nrow=total.return, ncol=FB)
     beta[1,] <- init.beta
   }else{
     beta <- matrix(init.beta, nrow=1, ncol=FB)
   }
   
   if ("bi"%in% monitor){ 
     bbb <- matrix(NA, nrow=total.return, ncol=RB*GU)
     bbb[1,] <- matrix(init.b, nrow=1, ncol=RB*GU)
   }else{
     bbb <- matrix(init.b, nrow=1, ncol=RB*GU)
   }

   if (timeint.add==1) {
     TPP <- TP+1
   }else{
     TPP <- TP
   }
   if ("ct"%in% monitor){ 
    ccc <- matrix(NA, nrow=total.return, ncol=GT*TPP)
    ccc[1,] <- matrix(init.c, nrow=1, ncol=GT*TPP)
   }else{
    ccc <- matrix(init.c, nrow=1, ncol=GT*TPP)
   }
   EDim <- (TPP^2-TPP)/2+TPP
   if ("E"%in% monitor){
     EEE <- matrix(NA, nrow=total.return, ncol=EDim)
     EEE[1,] <- matrix(0, nrow=1, ncol=EDim)
   }else{
     EEE <- matrix(0, nrow=1, ncol=EDim)
   }
   
   DDim <- ((RB+TP)^2-(RB+TP))/2+RB+TP
   if ("D"%in% monitor){
     DDD <- matrix(NA, nrow=total.return, ncol=DDim)
     DDD[1,] <- matrix(0, nrow=1, ncol=DDim)
   }else{
     DDD <- matrix(0, nrow=1, ncol=DDim)
   }

  if ("theta"%in% monitor){
     theta <- matrix(NA, nrow=total.return, ncol=TPP*GT)
     theta[1,] <- matrix(init.theta, nrow=1, ncol=TPP*GT)
   }else{
     theta <- matrix(init.theta, nrow=1, ncol=TPP*GT)
   }

   if ("gamma"%in% monitor){
     gamma <- matrix(NA, nrow=total.return, ncol=TP*GU)
     gamma[1,] <- matrix(init.gamma, nrow=1, ncol=TP*GU)
   }else{
     gamma <- matrix(init.gamma, nrow=1, ncol=TP*GU)
   }
     
   if ("u"%in% monitor){
     u.aux <- matrix(NA, nrow=total.return, ncol=N)
     u.aux[1,] <- matrix(rnorm(N), nrow=1, ncol=N)
   }else{
     u.aux <- matrix(rnorm(N), nrow=1, ncol=N)
   }
   
   if ("ystar"%in% monitor){ 
     ystarm <- matrix(NA, nrow=total.return, ncol=N)
     ystarm[1,] <- matrix(0, nrow=1, ncol=N)
   }else{
     ystarm <- matrix(0, nrow=1, ncol=N)
   }

  ##################################################
  # Calculate those quantities which will
  # not be updated in MCMC iterations
  ################################################## 
     invB <- pseudoinverse(B0)
     PBbeta <- invB%*%beta0
     En <- (e0+GT)/2
     En2 <- En*2
     if(is.matrix(E0)){
       invE <- pseudoinverse(E0)
     }
     Dn <- (e0+GT)/2
     Dn2 <- Dn*2
     if(is.matrix(D0)){
       invD <- pseudoinverse(D0)
     }
     cat("Data are arranged. MCMC begins.\n")
   
  ################################################
  ## MCMC updating loop begins here
  ################################################
     
  #Tracking the process of MCMC simulation
    mcmc <- (m+burnin)*thin
    for(g in 2:(mcmc)) {
      if (tracking>0){
       if ((g%% tracking ==0)| g==(m*thin)){
          cat("g=",g,"\n")
        }
     }
             
   # The current values of parameters

   if (g==2){
     current.bi <- matrix(bbb[1,], ncol=RB)           # current values of ci and beta3i
     current.theta <- matrix(theta[1,], ncol=TPP)     # current values of thetat
     current.gamma <- matrix(gamma[1,], ncol=TP)      # current values of gammai
     current.beta <- beta[1,]                         # current values of beta
     current.ystar <- ystarm[1,]
     current.phi <- init.phi
   }else{
     current.bi <- new.bi                            # current values of ci and beta3i
     currect.theta <- new.theta                      # current values of theta.t
     current.gamma <- new.gamma                      # current values of gamma.i
     current.beta <- new.beta                        # current values of beta
     current.ystar <- new.ystar                      # the current value of ystar
     current.phi <- new.phi
   }
  
  ###################################################
  # Block 1:  Covariance Matrices (random errors)
  ###################################################  

   if (TPP==1){ # if S only contains 1 element (most likely a constant)
     Ei<-rgamma(1,En, rate=(E0+sum(current.theta^2))/2)
     Ei <- as.matrix(Ei)
   }else{
     OEin <- invE+sumtrans(Mat=current.theta)
     sE<-rWishart (1, En2, pseudoinverse(OEin))
     Ei <- Wishart(element=sE, col=TPP)
   }

       
   current.unit <- cbind(current.bi, current.gamma)

   if (ncol(current.unit)==1){ # if W only contains 1 element (most likely a constant)
     Di<-rgamma(1,Dn, rate=(D0+sum(current.unit^2))/2)
     Di <- as.matrix(Di)
   }else{
     ODin <- invD+sumtrans(Mat=current.unit)
     sD<-rWishart (1, Dn2, pseudoinverse(ODin))
     Di <- Wishart(element=sD, col=(RB+TP))
  }

  #################################
  ### Prepare for data augmentation
  #################################
  # Compute the covariance-matrix of the serially correlated error \xi
    if (nlag==1){
      var.phig <- 1/(1- current.phi^2)                          
      TCov <- Cov.AR1(col=GT, rho= current.phi,sigma=var.phig)
    }else{
      CCov <- Cov.ARp(phig=current.phi, col=GT)        
      TCov <- CCov[[2]]
      Sig <- TCov[1:nlag, 1:nlag]
      SSig<-solve(Sig) 
      var <- Sig[1,1]
    }

   # Variables for outputs from the loop 
       
    betaterm1 <- matrix(0, nrow=FB, ncol=FB)     # for posterior covariance of beta
    betaterm2 <- rep(0, FB)                      # for posterior mean of beta
    new.u<- rep(NA, N)
    kap <- rep(NA, GU)                           # Store the kappa for variance
    qu <- rep(NA,N)                              # Store the output of the auxiliary term
    new.ystar <- rep(NA, N)
    Colum <- (RB+TP)
    new.unit2 <- matrix(NA, nrow=GU, ncol=Colum)
       
    ###################################################
    # Block 2:  data augmentation for Y and U (and prepare for beta )
    ###################################################
    for (n in 1: GU){                                 # loop for each unit
      rows <- entry[which(id1==n)]                    # the obserbations of unit n
      nT <- length(rows)                              # length of those observations
      yi <- Y[rows]                                   # Observed responses of n
      ystari <- current.ystar[rows]                   # current latent responses of n
      Xi <- as.matrix(X[rows,])
      Wi <- as.matrix(W[rows,])
      Si <- current.gamma[n,]
      if (timeint.add==1){
        Si <- c(1, Si)
      }
      bbbn <- current.bi[n,]
      TT <- id2[rows]
      LT <- max(TT)-min(TT)+1
      KIU <- min(TT):max(TT)
      AS <- KIU%in% TT

      Tccin <- as.matrix(current.theta[TT,])
      rant <- c()
      for (l in 1:nT){
        rant[l] <- Si%*%as.matrix(Tccin[l,])
      }

   #####################
   ## Updata U
   #####################
      nnCov <- as.matrix(TCov[1:LT, 1:LT])
      nCov <- as.matrix(nnCov[AS, AS])
      nVK <- V.matrix2(Omega=nCov, col=nT)
      nV <- nVK[[2]]                          
      nK <- nVK[[1]]                  
      Cov.U <- pseudoinverse(diag(nT) + 1/nK*t(nV)%*%nV)   
      mean.U <- Cov.U %*% t(nV) %*% (ystari-Xi%*%cbind(current.beta)-Wi%*%cbind(bbbn)-rant)/nK
      u.draw <- mvrnorm(1, mu=mean.U, Sigma=Cov.U)
      new.u[rows] <- u.draw
      u.decom <- nV%*%cbind(u.draw)
      qu[rows] <- u.decom
      kap[n] <- nK
      
  # Updata ystar
    mean.ystar <- Xi%*% cbind(current.beta)+Wi%*%cbind(bbbn)+rant+u.decom
    ystar <- sapply(c(1:nT),
             function(r){BTnormlog(x=yi[r], mu=mean.ystar[r], sigma=sqrt(nK))})
 
    new.ystar[rows] <- ystar

   # End of updating Y_i and U_i
   #*************************************************************#
   # 2) Begin to calculate the quantities for updating {b_i} and gamma
      
   if (n==1){
      if (timeint.add==1) Tccin <- as.matrix(Tccin[,-1])
      NW <- cbind(Wi, Tccin)
      Cov.b <- pseudoinverse(Di + t(NW)%*%pseudoinverse(nCov)%*%NW)
      mean.b <- Cov.b%*% t(NW)%*%pseudoinverse(nCov)%*%(ystar-Xi%*% cbind(current.beta))
      u1 <- mean.b[1:RB]
      sigma11 <- Cov.b[1:RB, 1:RB]
      x2 <- rep(1, TP)
      if (length(u1)==1){
        bb.draw <- rnorm(1, u1, sqrt(sigma11))
        b.drawn <- as.matrix(c(bb.draw, x2))
      }else{
        bb.drawn <- mvrnorm(1, mu=u1, Sigma=sigma11)
        b.drawn <- as.matrix(c(bb.drawn, x2))
      } 
    }else{
      if (timeint.add==1) Tccin <- as.matrix(Tccin[,-1])
        NW <- cbind(Wi, Tccin)
        Cov.b <- pseudoinverse(Di + t(NW)%*%pseudoinverse(nCov)%*%NW)
        mean.b <- Cov.b%*% t(NW)%*%pseudoinverse(nCov)%*%(ystar-Xi%*% cbind(current.beta))
      if (length(mean.b)==1){
        b.drawn <- as.matrix(rnorm(1, mean.b, sqrt(Cov.b)))
      }else{
          b.drawn <- mvrnorm(1, mu=mean.b, Sigma=Cov.b)
      }
    }
      new.unit2[n,] <- as.vector(b.drawn)
      
      nbbbn <- as.matrix(b.drawn[-c((Colum-TP+1):Colum)])
      ngamma <- as.matrix(b.drawn[c((Colum-TP+1):Colum)])
      nrant <- c()
      for (l in 1:nT){                                # Using loop to do it
        nrant[l] <- Tccin[l,]%*%ngamma
      }
           
   # calculate BB1 and BB2 for beta
      var.beta <- pseudoinverse(nCov)
      betaterm1 <- betaterm1+t(Xi) %*% var.beta %*%Xi
      betaterm2 <- betaterm2+t(Xi)%*% var.beta %*% (ystar-Wi%*% cbind(nbbbn)-nrant)
    }

       
     ############################
     ##Block 3: Beta and {b_i, \gamma_i}
     ############################
      # 1) {b_i, \gamma_i}
      new.bi <- as.matrix(new.unit2[, -c((Colum-TP+1):Colum)])
      new.gamma <- as.matrix(new.unit2[,  c((Colum-TP+1):Colum)])
       # 2) \beta
      B1<-pseudoinverse(betaterm1+invB)
      beta1<-t(B1%*%(betaterm2+PBbeta))
      new.beta<-mvrnorm(1, beta1, B1)
      # End of Block 3
       
    #######################################
    # Block 4: Time-Specific Random Effect
    ######################################
      new.theta <-matrix(NA, nrow=GT, ncol=TPP)
      for (j in 1: GT){
        row1 <- entry[which(id2==j)]    # Those entries with t=j
        H <- length(row1)               #  #observations at t period
        iden1 <- id1[row1]              # The units at t period
        tttn <- as.matrix(new.bi[iden1,])  # The random unit-specific effects of
        yt1 <- new.ystar[row1]
        Xt1 <- as.matrix(X[row1,])
        Wt1<- as.matrix(W[row1,])
        St1 <- as.matrix(new.gamma[iden1,])
        if(timeint.add==1){
          St1 <- cbind(1, St1)
        }
        Ut1 <- qu[row1]
        kapt1 <- kap[iden1]
        
       term7 <- c()
       for (q in 1:H){
         term7[q] <- Wt1[q,]%*% cbind(tttn[q,])
       }
            
       PM <- matrix(0, nrow=H, ncol=H)
       for (c in 1:H){
         PM[c,c] <- kapt1[c]
       }
            
       cov.S <- pseudoinverse(Ei+t(St1)%*%pseudoinverse(PM)%*%St1)
       mean.S <- cov.S%*% t(St1) %*%pseudoinverse(PM)%*%(yt1-term7-Xt1%*%cbind(new.beta)-Ut1)
       if (length(mean.S)==1){
         supd <- as.matrix(rnorm(1, mean.S, sqrt(cov.S)))
       }else{
         supd <- mvrnorm(1, mean.S, cov.S)
       }
       new.theta[j,] <- as.vector(supd)
      }
 
  # End of Block 4
       
    ###############################
    ## Block 5: Autoregressive Coef.s
    ###############################
      time.mutiple <- c()
      for (n in 1: GU){                                # loop for each unit
        St1 <- new.gamma[n,]
        if (timeint.add==1) St1 <- c(1, St1)
        rows <- entry[which(id1==n)]                    # the obserbations of unit n
        TT <- id2[rows]
        Tccin <- as.matrix(new.theta[TT,])
        time.mutiple[rows] <-  St1%*%t(Tccin)
      }
      time.mutiple <- as.matrix(time.mutiple)

      if(any(id1!=id3)){
        new.bii <- matrix(NA, nrow=GUU, ncol=ncol(new.bi))
        new.bii[1,] <- new.bi[1,]
        for (h in 2: GUU){
         rowid3 <- entry[id3==h]
         hh <- unique(id1[rowid3])
         new.bii[h,] <- new.bi[hh,]
        }
      }else{
        new.bii <- new.bi
      }
      ###### MH Algorithm for AR(1) ########
      if (nlag==1){
       PhiPr <- PhiTerm0(Entry=entry, unit=id3, NN=GU, yystar=new.ystar,newb=new.bii, PX=X,
                         PW=W, PS=time.mutiple, Ti=id2, CC=as.matrix(rep(1, nrow(time.mutiple))),
                         newbeta=new.beta)
       PhiSig <- PhiPr[1]
       PhiMean <- PhiPr[2]
       hatPhi<-1/PhiSig
       hatphi<-hatPhi*PhiMean
       phi.draw <- rtnorm(1, hatphi, sqrt(hatPhi), lower=-1, upper=1)
       if (tracking>0){
         if ((g%% tracking ==0)| g==(m*thin)){       
          cat("phi.draw ", phi.draw, "\n", sep=" ")
         }
       }
       # Accept or Reject
       
       alpha<-AcceptRate0(rho=phi.draw, orho=current.phi, unit=id3, precision=1,AY=new.ystar,
                          AX=X, AW=W, AS=time.mutiple, Abeta=new.beta, Ab=new.bii, AT=id2,
                          Ac=as.matrix(rep(1, nrow(time.mutiple))), Entry=entry)

        u <- runif(1)
        if( u <= alpha ) {
            new.phi<-phi.draw
        }else{
           new.phi<-current.phi
        }
     }

      ###### MH Algorithm for AR(p>1) ########      
      if (nlag>1){ 
         PhiPr <- PhiTermARpT(Entry=entry, unit=id3, NN=GU, yystar=new.ystar, newb=new.bii, PX=X,
                              PW=W, PS=time.mutiple, Ti=id2, CC=as.matrix(rep(1, nrow(time.mutiple))),
                              newbeta=new.beta, nlag=nlag)
         PhiSig <- PhiPr[,-(nlag+1)]
         PhiMean <- PhiPr[,(nlag+1)]
         phigg<-PhiPropARp(Ymt=PhiSig, Ymt2=PhiMean,scale=1)
     if (tracking>0){
      if ((g%% tracking ==0)| g==(m*thin)){         
        cat("phi.draw ", phigg, "\n", sep=" ")
      }
     }    
     # Accept or Reject 
         
        alpha<-AcceptRateARpT(nrho=phigg, orho=current.phi, Osig=Sig , unit=id3, AY=new.ystar,
                              AX=X, AW=W, AS=time.mutiple, Abeta=new.beta, Ab=new.bii,  AT=id2,
                              Ac=as.matrix(rep(1, nrow(time.mutiple))), Entry=entry, f=g, nlag=nlag)
     #   cat("alpha ", alpha, "\n", sep=" ")
       u <- runif(1)
        if( u <= alpha ) {
            new.phi <- phigg
        } else {
            new.phi<- current.phi
        }
     }

   ############################################
   #Store the updated values
   ############################################       
  if ((g%%thin ==0)| g==m){
         gg <- g%/%thin
         if ("ystar"%in% monitor){
           ystarm[gg,] <- new.ystar
         }
         
         if ("u"%in% monitor){
           u.aux[gg,] <- new.u
         }
       
         if ("bi" %in% monitor){
           bbb[gg,] <- as.vector(new.bi)
         }
        
         
         if ("gamma" %in% monitor){
           gamma[gg,] <- as.vector(new.gamma)
         }
        
         if ("theta"%in% monitor){
           theta2 <- as.vector(new.theta)
           theta[gg,] <- theta2
         }
         
         if ("beta"%in% monitor){
           beta[gg,] <- new.beta
         }
       
         if ("rho"%in% monitor){
           if (nlag==1){
             phi[gg] <- new.phi
           }else{
             phi[gg,] <- new.phi
           }
         }
       
         if ("D"%in% monitor){
           if ((RB+TP)==1){
             DDD[gg,] <-  as.vector(Di)
           }else{
             DDD[gg,] <-  sD
           }
          }
       
         if ("E"%in% monitor){
           if (TPP==1){
           EEE[gg,] <- as.vector(Ei)
         }else{
           EEE[gg,] <- sE
           }
         }
       }
    }

     ###############################
     # Return the MCMC results
     ###############################

     if ("rho"%in% monitor){
       if (nlag==1){posterior.phi <- phi[-c(1:burnin)]
        }else{
          posterior.phi <- phi[-c(1:burnin),]
        }
       }else{
        posterior.phi <- phi
      }
     if ("beta"%in% monitor){ 
       posterior.beta <- beta[-c(1:burnin),]
     }else{
       posterior.beta <- beta
     }
     if ("bi"%in% monitor){ 
       posterior.unit.res <- bbb[-c(1:burnin),]
     }else{
       posterior.unit.res <- bbb
     }
     if ("gamma"%in% monitor){ 
       posterior.timeunit.res <- gamma[-c(1:burnin),]
     }else{
       posterior.timeunit.res <- gamma
     }
     if ("theta"%in% monitor){ 
       posterior.timetime.res <- theta[-c(1:burnin),]
     }else{
       posterior.timetime.res <- theta
     }
     
    if ("E"%in% monitor){ 
      posterior.time.cov <- EEE[-c(1:burnin),]
    }else{
      posterior.time.cov <- EEE
    }
     
    if ("D"%in% monitor){ 
      posterior.unit.cov <- DDD[-c(1:burnin),]
    }else{
      posterior.unit.cov <- DDD
    }

     if ("u"%in% monitor){
      posterior.u.aux <- u.aux[-c(1:(burnin/thin)),]
     }else{
      posterior.u.aux <- u.aux
     }
     
     if ("ystar"%in% monitor){ 
       posterior.data <- ystarm[-c(1:(burnin/thin)),]
     }else{
       posterior.data <- ystarm
     }

    return(list(Phi=posterior.phi, Beta=posterior.beta,  bi=posterior.unit.res,
                gamma=posterior.timeunit.res, theta=posterior.timetime.res,
                var.theta=posterior.time.cov, Cov.betai=posterior.unit.cov,
                auxiliary.u=posterior.u.aux, data.augmented=posterior.data))
   }




#################################
## AR0
#################################

  MulFac.AR0 <- function(y,X1, Wi, Ai, Ut, TP, Unit.Index, Time.Index,  m=10000, burnin=5000,
                        init.ystar=0, init.beta=0,unitint.add=0,timeint.add=0, init.b=0, init.c,
                         init.theta=0, init.gamma=0, beta0, B0, D0, d0, E0, e0,  nlag=0,
                          monitor=c( "beta", "bi", "ct", "D", "E", "ystar", "u", "theta", "gamma"),
                          thin=1, tracking){
    ############################
     #y \in {0, 1} error checking
     ############################
     yy <- y[!is.na(y)]
     if (sum(yy!=0 & yy!=1) > 0) {
       cat("Elements of y equal to something other than 0 or 1.\n")
       stop("Check data and call MulFac.AR() again. \n") 
     }
     
    #****************** Arrange Data for MCMC updating  **********************#
     
     # To use the RegressionData for this model specification, the arguments of
     #  Ft and St are set to be certain values since the data inpute does not contain
     #  information for Ft and St
     Ft <- "NULL"
     St <- as.matrix(rep(1, length(y)))
     # Call RegressionData() and rearrange the data to fit the reduced model
     NewData <- RegressionData(yob=y, Unit=Unit.Index, Time=Time.Index, Fixed=X1,
                               UnitRandom=Wi, TimeRandom=St,
                               UnitPred=Ai, TimePred=Ut)

    Y <- NewData[[3]]
    entry1 <- 1: length(Y)
    newid1 <- NewData[[1]]
    newid2 <- NewData[[2]]
    id11 <- Index(newid1)     #the unit diex
    uniqunit <- as.numeric(table(id11))   # how many time periods for each unit
    jj <- which(uniqunit<= nlag)
    if (length(jj)>0){						
      entry2 <- entry1[id11==jj]
      Y <- Y[-entry2]
      X <- as.matrix(NewData[[4]][-entry2,])
      WP <- as.matrix(NewData[[5]][-entry2,])
    }else{
      X <- as.matrix(NewData[[4]])
      WP <- as.matrix(NewData[[5]])
    }
    WPP <- ncol(WP)
    KW <- ncol(Ut)
    FB <- ncol(X)
    W <- as.matrix(cbind(WP, X[, (FB-KW+1):FB]))
    if (unitint.add==TRUE){
      W <- cbind(1, W)
    }
   N <- length(Y)
     
   # Check the dimension of the parameter vectors
     
   #beta: fixed parameter vector 
   PFB <- length(beta0)
    if (FB!=PFB){
       cat("Length of the prior on fixed effects is incorrect.\n")
       cat("Length of fixed effect should be ", FB, "\n", sep=" ")
       stop("Respecify beta0 and B0, call GLMMARp.Binary() again. \n") 
    }
     
   # bi: random fixed parameter vector
     RB <- ncol(W)
     if(is.matrix(D0)){ 
       PRB <- ncol(D0)-TP
     }else{
       PRB <- length(D0)-TP
     }
   if (RB!=PRB){
       cat("Dimension of the prior on unit random effects is incorrect.\n")
       cat("Length of D0 should be ", RB+TP, "*", RB+TP,  "\n", sep=" ")
       stop("Respecify D0, call GLMMARp.Binary() again. \n") 
    }

  # in case of unbalanced TSCS data structure,
   # create new indices to accomendate such a situation 						
   entry <- seq(1: N)                   # the global length of observations
   if (length(jj)>0){
     # unit id repeating t times t=# observations
     id1 <- Index(id11[-entry2])                
     # year id repeating n times n=# unit in each time period
     id2 <- (newid2[-entry2])-min(newid2[-entry2])+1        
     GU <- max(id1)                       # how many units
     GT <- max(id2)                       # how many time periods in total
     id3 <- id1
     GU <- length(unique(id3))
     n=1
     while (n <= length(unique(id3))){
       rows <- entry[which(id3==n)]                   
       nT <- length(rows)
       TT <- id2[rows]
       if (nT!=1){
         for (j in 2:nT){
          GHK <- TT[j]-TT[j-1]
          if (GHK >1){
            PL <- rows[j]
            id3[PL:N] <- id3[PL:N]+1
          }
        }
       }
       n <- n+1
     }

    entry3 <- 1: length(Y)
    uniqunit <- as.numeric(table(id3))
    jj <- which(uniqunit<= nlag)
    entry4 <- entry1[id3==jj]
    if(length(entry4)>0){
      Y <- Y[-entry4]
      X <- as.matrix(X[-entry4,])
      W <- as.matrix(W[-entry4,])
      id1 <- Index(id1[-entry4])
      id2 <- Index(id2[-entry4])
      id3 <- Index(id3[-entry4])
    }
   }else{
     id1 <- id11
     id2 <- Index(newid2)
   }
    id1 <- Index(id1)                 # unit id repeating t times t=# observations
    id2 <- (id2)-min(id2)+1        # year id repeating n times n=# unit in each time period
    id3 <- id1						
    GU <- max(id1)                       # how many units
    GT <- max(id2)
    GUU <- length(unique(id3))

							  
     #############################
     # Storage and Initial Values
     #############################
     total.return <- m+burnin

     if ("beta"%in% monitor){ 
       beta <- matrix(NA, nrow=total.return, ncol=FB)
       beta[1,] <- init.beta
     }else{
       beta <- matrix(init.beta, nrow=1, ncol=FB)
     }
   
     if ("bi"%in% monitor){ 
       bbb <- matrix(NA, nrow=total.return, ncol=RB*GU)
       bbb[1,] <- matrix(init.b, nrow=1, ncol=RB*GU)
     }else{
       bbb <- matrix(init.b, nrow=1, ncol=RB*GU)
     }

     if (timeint.add==1) {
       TPP <- TP+1
     }else{
       TPP <- TP
     }
     
     if ("ct"%in% monitor){ 
      ccc <- matrix(NA, nrow=total.return, ncol=GT*TPP)
      ccc[1,] <- matrix(init.c, nrow=1, ncol=GT*TPP)
     }else{
      ccc <- matrix(init.c, nrow=1, ncol=GT*TPP)
     }
     
     EDim <- (TPP^2-TPP)/2+TPP
     if ("E"%in% monitor){
      EEE <- matrix(NA, nrow=total.return, ncol=EDim)
      EEE[1,] <- matrix(0, nrow=1, ncol=EDim)
     }else{
      EEE <- matrix(0, nrow=1, ncol=EDim)
     }
   
     DDim <- ((RB+TP)^2-(RB+TP))/2+RB+TP
     if ("D"%in% monitor){
       DDD <- matrix(NA, nrow=total.return, ncol=DDim)
       DDD[1,] <- matrix(0, nrow=1, ncol=DDim)
     }else{
       DDD <- matrix(0, nrow=1, ncol=DDim)
     }

    if ("theta"%in% monitor){
      theta <- matrix(NA, nrow=total.return, ncol=TPP*GT)
      theta[1,] <- matrix(init.theta, nrow=1, ncol=TPP*GT)
    }else{
     theta <- matrix(init.theta, nrow=1, ncol=TPP*GT)
    }

   if ("gamma"%in% monitor){
     gamma <- matrix(NA, nrow=total.return, ncol=TP*GU)
     gamma[1,] <- matrix(init.gamma, nrow=1, ncol=TP*GU)
   }else{
     gamma <- matrix(init.gamma, nrow=1, ncol=TP*GU)
   }
   
   if ("ystar"%in% monitor){ 
     ystarm <- matrix(NA, nrow=total.return, ncol=N)
     ystarm[1,] <- matrix(0, nrow=1, ncol=N)
   }else{
     ystarm <- matrix(0, nrow=1, ncol=N)
   }
							  
  ##################################################
  # Calculate those quantities which will
  # not be updated in MCMC iterations
  ################################################## 
     invB <- pseudoinverse(B0)
     PBbeta <- invB%*%beta0
     En <- (e0+GT)/2
     En2 <- En*2
     if(is.matrix(E0)){
       invE <- pseudoinverse(E0)
     }
     Dn <- (e0+GT)/2
     Dn2 <- Dn*2
     if(is.matrix(D0)){
       invD <- pseudoinverse(D0)
     }
     cat("Data are arranged. MCMC begins.\n")
     
  ################################################
  ## MCMC updating loop begins here
  ################################################
     
  #Tracking the process of MCMC simulation
    mcmc <- (m+burnin)*thin
    for(g in 2:(mcmc)) {
      if (tracking>0){
       if ((g%% tracking ==0)| g==(m*thin)){
          cat("g=",g,"\n")
        }
     }
             
   # The current values of parameters

   if (g==2){
     current.bi <- matrix(bbb[1,], ncol=RB)          # current values of ci and beta3i
     current.theta <- matrix(theta[1,], ncol=TPP)      # current values of thetat
     current.gamma <- matrix(gamma[1,], ncol=TP)      # current values of gammai
     current.beta <- beta[1,]                         # current values of beta
     current.ystar <- ystarm[1,]
   }else{
     current.bi <- new.bi                            # current values of ci and beta3i
     currect.theta <- new.theta                      # current values of theta.t
     current.gamma <- new.gamma                      # current values of gamma.i
     current.beta <- new.beta                        # current values of beta
     current.ystar <- new.ystar                      # the current value of ystar
   }
  
  ###################################################
  # Block 1:  Covariance Matrices (random errors)
  ###################################################  

   if (TPP==1){ # if S only contains 1 element (most likely a constant)
     Ei<-rgamma(1,En, rate=(E0+sum(current.theta^2))/2)
     Ei <- as.matrix(Ei)
   }else{
     OEin <- invE+sumtrans(Mat=current.theta)
     sE<-rWishart (1, En2, pseudoinverse(OEin))
     Ei <- Wishart(element=sE, col=TPP)
   }

       
   current.unit <- cbind(current.bi, current.gamma)

   if (ncol(current.unit)==1){ # if W only contains 1 element (most likely a constant)
     Di<-rgamma(1,Dn, rate=(D0+sum(current.unit^2))/2)
     Di <- as.matrix(Di)
   }else{
     ODin <- invD+sumtrans(Mat=current.unit)
     sD<-rWishart (1, Dn2, pseudoinverse(ODin))
     Di <- Wishart(element=sD, col=(RB+TP))
  }
  
    ###################################################
    # Block 2:  data augmentation for Y (and prepare for beta )
    ###################################################

    betaterm1 <- matrix(0, nrow=FB, ncol=FB)     # for posterior covariance of beta
    betaterm2 <- rep(0, FB)                      # for posterior mean of beta 
    new.ystar <- rep(NA, N)
    Colum <- (RB+TP)
    new.unit2 <- matrix(NA, nrow=GU, ncol=Colum)
       
    # Big loop for data augmentation and values for beta 
    for (n in 1: GU){                                 # loop for each unit
      rows <- entry[which(id1==n)]                    # the obserbations of unit n
      nT <- length(rows)                              # length of those observations
      yi <- Y[rows]                                   # Observed responses of n
      ystari <- current.ystar[rows]                   # current latent responses of n
      Xi <- as.matrix(X[rows,])
      Wi <- as.matrix(W[rows,])
      Si <- current.gamma[n,]
      if (timeint.add==1){
        Si <- c(1, Si)
      }
      bbbn <- current.bi[n,]
      TT <- id2[rows]
      LT <- max(TT)-min(TT)+1
      KIU <- min(TT):max(TT)
      AS <- KIU%in% TT

      Tccin <- as.matrix(current.theta[TT,])
      rant <- c()
      for (l in 1:nT){
        rant[l] <- Si%*%as.matrix(Tccin[l,])
      }

    # Updata ystar
    mean.ystar <- Xi%*% cbind(current.beta)+Wi%*%cbind(bbbn)+rant
    ystar <- sapply(c(1:nT),
             function(r){BTnormlog(x=yi[r], mu=mean.ystar[r], sigma=1)})
 
    new.ystar[rows] <- ystar


   #*************************************************************#
    #  Begin to calculate the quantities for updating {ct,f} and gamma
   if (n==1){
      if (timeint.add==1) Tccin <- as.matrix(Tccin[,-1])
      NW <- cbind(Wi, Tccin)
      Cov.b <- pseudoinverse(Di + t(NW)%*%NW)
      mean.b <- Cov.b%*% t(NW)%*%(ystar-Xi%*% cbind(current.beta))
      u1 <- mean.b[1:RB]
      sigma11 <- Cov.b[1:RB, 1:RB]
      x2 <- rep(1, TP)
      if (length(u1)==1){
        bb.draw <- rnorm(1, u1, sqrt(sigma11))
        b.drawn <- as.matrix(c(bb.draw, x2))
      }else{
        bb.drawn <- mvrnorm(1, mu=u1, Sigma=sigma11)
        b.drawn <- as.matrix(c(bb.drawn, x2))
      } 
    }else{
      if (timeint.add==1) Tccin <- as.matrix(Tccin[,-1])
      NW <- cbind(Wi, Tccin)
      Cov.b <- pseudoinverse(Di + t(NW)%*%NW)
      mean.b <- Cov.b%*% t(NW)%*%(ystar-Xi%*% cbind(current.beta))
      if (length(mean.b)==1){
        b.drawn <- as.matrix(rnorm(1, mean.b, sqrt(Cov.b)))
      }else{
        b.drawn <- mvrnorm(1, mu=mean.b, Sigma=Cov.b)
     }
   }
   new.unit2[n,] <- as.vector(b.drawn)
      
   nbbbn <- as.matrix(b.drawn[-c((Colum-TP+1):Colum)])
   ngamma <- as.matrix(b.drawn[c((Colum-TP+1):Colum)])
   nrant <- c()
   for (l in 1:nT){                               
     nrant[l] <- Tccin[l,]%*%ngamma
   }
           
   # calculate BB1 and BB2 for beta
    
    betaterm1 <- betaterm1+t(Xi) %*%Xi
    betaterm2 <- betaterm2+t(Xi)%*% (ystar-Wi%*% cbind(nbbbn)-nrant)
  }

       
     ###################################
     ##Block 3: Beta and {b_i, \gamma_i}
     ###################################
      # 1) {b_i, \gamma_i}
     new.bi <- as.matrix(new.unit2[, -c((Colum-TP+1):Colum)])
     new.gamma <- as.matrix(new.unit2[,  c((Colum-TP+1):Colum)])

     # 2) \beta
     B1<-pseudoinverse(betaterm1+invB)
     beta1<-t(B1%*%(betaterm2+PBbeta))
     new.beta<-mvrnorm(1, beta1, B1)
       
    #######################################
    # Block 4: Time-Specific Random Effect
    ######################################
       
    new.theta <-matrix(NA, nrow=GT, ncol=TPP)
    for (j in 1: GT){
      row1 <- entry[which(id2==j)]    # Those entries with t=j
      H <- length(row1)               #  #observations at t period
      iden1 <- id1[row1]              # The units at t period
      tttn <- as.matrix(new.bi[iden1,])         # The random unit-specific effects of
      yt1 <- new.ystar[row1]
      Xt1 <- as.matrix(X[row1,])
      Wt1<- as.matrix(W[row1,])
      St1 <- as.matrix(new.gamma[iden1,])
      if(timeint.add==1){
        St1 <- cbind(1, St1)
      }
            
      term7 <- c()
      for (q in 1:H){
        term7[q] <- Wt1[q,]%*% cbind(tttn[q,])
      }
      cov.S <- pseudoinverse(Ei+t(St1)%*%St1)
      mean.S <- cov.S%*% t(St1)%*%(yt1-term7-Xt1%*%cbind(new.beta))
      if (length(mean.S)==1){
       supd <- as.matrix(rnorm(1, mean.S, sqrt(cov.S)))
      }else{
       supd <- rmvnorm(1, mean.S, cov.S)
      }
      new.theta[j,] <- as.vector(supd)
    }
 
  # End of Block 4
  
   ############################################
   #Store the updated values
   ############################################        
  if ((g%%thin ==0)| g==m){
         gg <- g%/%thin
         if ("ystar"%in% monitor){
           ystarm[gg,] <- new.ystar
         }

       
         if ("bi" %in% monitor){
           bbb[gg,] <- as.vector(new.bi)
         }
        
         
         if ("gamma" %in% monitor){
           gamma[gg,] <- as.vector(new.gamma)
         }
        
         if ("theta"%in% monitor){
           theta2 <- as.vector(new.theta)
           theta[gg,] <- theta2
         }
         
         if ("beta"%in% monitor){
           beta[gg,] <- new.beta
         }
       
         if ("D"%in% monitor){
           if ((RB+TP)==1){
             DDD[gg,] <-  as.vector(Di)
           }else{
             DDD[gg,] <-  sD
           }
          }
       
         if ("E"%in% monitor){
           if (TPP==1){
           EEE[gg,] <- as.vector(Ei)
         }else{
           EEE[gg,] <- sE
           }
         }
       }
    }

     ###############################
     # Return the MCMC results
     ###############################  
   if ("beta"%in% monitor){ 
     posterior.beta <- beta[-c(1:burnin),]
   }else{
     posterior.beta <- beta
   }
   if ("bi"%in% monitor){ 
     posterior.unit.res <- bbb[-c(1:burnin),]
   }else{
     posterior.unit.res <- bbb
   }
   if ("gamma"%in% monitor){ 
     posterior.timeunit.res <- gamma[-c(1:burnin),]
   }else{
     posterior.timeunit.res <- gamma
   }
   if ("theta"%in% monitor){ 
     posterior.timetime.res <- theta[-c(1:burnin),]
   }else{
     posterior.timetime.res <- theta
   } 
   if ("E"%in% monitor){ 
     posterior.time.cov <- EEE[-c(1:burnin),]
   }else{
     posterior.time.cov <- EEE
   }    
   if ("D"%in% monitor){ 
     posterior.unit.cov <- DDD[-c(1:burnin),]
   }else{
     posterior.unit.cov <- DDD
   }
   if ("ystar"%in% monitor){ 
     posterior.data <- ystarm[-c(1:(burnin/thin)),]
   }else{
     posterior.data <- ystarm
   }

    return(list(Beta=posterior.beta,  bi=posterior.unit.res, gamma=posterior.timeunit.res,
                theta=posterior.timetime.res, var.theta=posterior.time.cov,
                Cov.betai=posterior.unit.cov, data.augmented=posterior.data))
   }



